SCBI Data
# From forestecology package
# devtools::install_github("rvalavi/blockCV")
library(blockCV)
library(sfheaders)
scbi_2013 <- read.csv("data/scbi.stem2.csv") %>%
select(treeID, stemID, sp, quadrat, gx, gy, dbh, date = ExactDate, codes, status) %>%
mutate(date = lubridate::mdy(date)) %>%
filter(gx < 300, gy > 300, gy < 600)
scbi_2018 <- read.csv("data/scbi.stem3.csv") %>%
select(treeID, stemID, sp, quadrat, gx, gy, dbh, date = ExactDate, codes, status) %>%
mutate(date = lubridate::mdy(date),
dbh = as.numeric(dbh)) %>%
filter(gx < 300, gy > 300, gy < 600)
census_df1 <- scbi_2013
census_df2 <- scbi_2018
id <- "stemID"
cluster_names <- data.frame(group = c(1:7), group_name = c("high_sla", "low_sla", "shrubs", "tall_light_wood", "evergreens", "oaks", "other"), group_alias = c("High SLA trees", "Low SLA trees", "Shrubs", "Tall trees - Light wood", "Evergreens", "Oaks", "Other"))
spptable <- read_csv("data/scbi.spptable.csv") %>%
left_join(cluster_names, by = "group") %>%
mutate(group = ifelse(is.na(group), 7, group),
group_name = ifelse(group == 7, "other", group_name),
group_alias = ifelse(group == 7, "Other", group_alias))
scbi_growth_df <-
# Merge both censuses and compute growth:
compute_growth(census_df1, census_df2, id) %>%
# Convert growth from cm to mm to make result comparable
mutate(growth = growth/10,
sp = as.factor(sp)) %>%
inner_join(spptable, by = "sp") %>%
filter(status == "A")# Add spatial information
cv_fold_size <- 100
max_dist <- 7.5
scbi_study_region <-
#tibble(x = c(0,400,400,0,0), y = c(0,0,640,640,0)) %>%
tibble(x = c(0,300,300,0,0), y = c(300,300,600,600,3)) %>%
sfc_polygon()
# Add buffer variable to data frame
scbi_growth_df <- scbi_growth_df %>%
add_buffer_variable(direction = "in", size = max_dist, region = scbi_study_region)
scbi_cv_grid <- spatialBlock(
speciesData = scbi_growth_df, theRange = 100, yOffset = 0.9999, k = 9, verbose = FALSE)
# Add foldID to data
scbi_growth_df <- scbi_growth_df %>%
mutate(
foldID = scbi_cv_grid$foldID
)
# Visualize grid
scbi_cv_grid$plots +
geom_sf(data = scbi_growth_df, aes(col=factor(foldID)), size = 0.1)
scbi_cv_grid_sf <- scbi_cv_grid$blocks %>%
st_as_sf()Create focal vs. competitor data
focal_vs_comp_scbi <- scbi_growth_df %>%
create_focal_vs_comp(max_dist, cv_grid_sf = scbi_cv_grid_sf, id = "stemID")
focal_vs_comp_scbi <- focal_vs_comp_scbi %>%
left_join(select(spptable, sp, focal_group = group_name), by = c("focal_sp" = "sp")) %>%
left_join(select(spptable, sp, comp_group = group_name), by = c("comp_sp" = "sp"))Spread into wide form for modeling (note: should we do this by genus?)
focal_vs_comp_scbi_wide <- focal_vs_comp_scbi %>%
group_by(focal_ID, focal_group, dbh, foldID, growth, comp_group) %>%
summarize(comp_basal_area = sum(comp_basal_area)) %>%
pivot_wider(names_from = comp_group, values_from = comp_basal_area, values_fill = 0)
group_order <- focal_vs_comp_scbi_wide %>%
ungroup() %>%
filter(focal_group != "other") %>%
count(focal_group, sort = T)
# set factor order by largest group
focal_vs_comp_scbi_wide <- focal_vs_comp_scbi_wide %>%
mutate(focal_group = factor(focal_group, levels = c(as.character(group_order$focal_group), "other")))
head(focal_vs_comp_scbi_wide)## # A tibble: 6 x 12
## # Groups: focal_ID, focal_group, dbh, foldID, growth [6]
## focal_ID focal_group dbh foldID growth high_sla low_sla shrubs
## <int> <fct> <dbl> <int> <dbl> <dbl> <dbl> <dbl>
## 1 4 low_sla 136. 6 0.103 2.34 4.58 1.96
## 2 5 shrubs 88 6 0.150 1.47 9.88 6.89
## 3 79 tall_light… 477. 9 -0.161 13.8 34.8 4.52
## 4 80 shrubs 51.5 6 0.253 0 0.356 2.45
## 5 96 other 23 9 0.262 5.09 0 1.53
## 6 101 tall_light… 654. 7 0.552 2.83 0.859 0
## # … with 4 more variables: tall_light_wood <dbl>, other <dbl>, oaks <dbl>,
## # evergreens <dbl># write_csv(focal_vs_comp_scbi_wide, "data/focal_vs_comp_scbi_wide.csv")Model
model <- lm(growth ~ dbh + focal_group + high_sla + low_sla + tall_light_wood + evergreens + oaks + shrubs + other, data = focal_vs_comp_scbi_wide)
summary(model)##
## Call:
## lm(formula = growth ~ dbh + focal_group + high_sla + low_sla +
## tall_light_wood + evergreens + oaks + shrubs + other, data = focal_vs_comp_scbi_wide)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.57235 -0.06850 -0.01825 0.04713 2.34239
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.146e-02 6.274e-03 12.984 < 2e-16 ***
## dbh 5.346e-04 1.295e-05 41.292 < 2e-16 ***
## focal_grouphigh_sla -2.108e-02 5.922e-03 -3.559 0.000375 ***
## focal_grouplow_sla -3.323e-02 6.127e-03 -5.424 6.05e-08 ***
## focal_groupshrubs -2.605e-03 6.974e-03 -0.374 0.708766
## focal_groupoaks -6.129e-02 8.991e-03 -6.816 1.02e-11 ***
## focal_groupevergreens 1.106e-01 1.802e-02 6.134 9.07e-10 ***
## focal_groupother -5.637e-03 6.095e-03 -0.925 0.355133
## high_sla -1.247e-05 1.466e-04 -0.085 0.932230
## low_sla -6.453e-05 1.060e-04 -0.609 0.542575
## tall_light_wood -3.430e-04 7.381e-05 -4.646 3.45e-06 ***
## evergreens -8.056e-04 4.491e-04 -1.794 0.072878 .
## oaks -5.625e-04 8.943e-05 -6.290 3.39e-10 ***
## shrubs 7.497e-03 1.447e-03 5.180 2.29e-07 ***
## other -7.219e-04 4.913e-04 -1.469 0.141836
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1413 on 6302 degrees of freedom
## Multiple R-squared: 0.3035, Adjusted R-squared: 0.3019
## F-statistic: 196.1 on 14 and 6302 DF, p-value: < 2.2e-16\[ y = \beta_0 + \beta_1d_f + \beta_2g_f + \sum_{g=1}^G \beta_{g}m_g + \varepsilon \\ \text{ where } d_f \text{ is the diameter at breast height of the focal tree, } g_f \text{ is the cluster group of the focal tree, } \\ g \in \{\text{tall trees/light wood, high SLA trees, low SLA trees, shrubs, oaks, evergreens, other}\}, \\ \text{ and } m_g \text{ is the biomass of all competitor trees of the specific focal group. Coefficients are with respect to tall trees/light wood} \]
Plot all trees in SCBI plot:
scbi_growth_subset <- scbi_growth_df %>%
filter(stemID %in% focal_vs_comp_scbi_wide$focal_ID)
ggplot(data = scbi_growth_subset) +
geom_sf(aes(col = sp), size = 0.5) +
labs(title = "All trees in SCBI site")
ggplot(data = scbi_growth_subset) +
geom_sf(aes(col = Genus), size = 0.5) +
labs(title = "All trees in SCBI site")
ggplot(data = scbi_growth_subset) +
geom_sf(aes(col = Family), size = 0.5) +
labs(title = "All trees in SCBI site")
ggplot(data = scbi_growth_subset) +
geom_sf(aes(col = factor(group_alias)), size = 0.5) +
labs(title = "All trees in SCBI site")
Plot all the trees, specific species
ggplot(data = scbi_growth_df %>% filter(sp == "nysy")) +
geom_sf(aes(col = sp), size = 0.5)+
labs(title = "All Nyssa sylvatica trees in Michigan Big Woods site")
ggplot(data = scbi_growth_df %>% filter(sp == "astr")) +
geom_sf(aes(col = sp), size = 0.5)+
labs(title = "All Asimina triloba trees in Michigan Big Woods site")
maples <- read_csv("https://rudeboybert.github.io/SDS390/static/data/maples.csv")Model loos like this: \[ \begin{eqnarray*} y &=& \beta_0 + \beta_{\text{dbh}}\cdot \text{dbh} + \beta_{\text{evergreen_bm}}\cdot \text{evergreen_bm} + \beta_{\text{maple_bm}}\cdot \text{maple_bm} + \beta_{\text{misc_bm}}\cdot \text{misc_bm}\\ && + \beta_{\text{oak}}\cdot \text{oak} + \beta_{\text{short_tree}}\cdot \text{short_tree} + \beta_{\text{shrub}}\cdot \text{shrub} + \epsilon \end{eqnarray*} \]
where
- \(\text{dbh}\) is the “starting” DBH in 2008
- \(\text{bm}\) is “biomass” as estimated by the basal area.